{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains various convergence tests for MPFA.\n",
    "\n",
    "They are included mainly to complement the tests for MPSA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.sparse.linalg import spsolve\n",
    "import scipy.sparse as sps\n",
    "import sympy\n",
    "from math import pi\n",
    "\n",
    "import setup_grids\n",
    "\n",
    "from porepy.numerics.fv import mpfa, fvutils\n",
    "\n",
    "from porepy.params.tensor import SecondOrderTensor as perm\n",
    "from porepy.params import bc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian errors - no perturbations\n",
      "[0.09744623673648027, 0.02563304567729676, 0.006512296044445567] pressure error\n",
      "[0.11792041665757692, 0.038488831450398144, 0.011371514390600216] flux error\n",
      "Cartesian errors - perturbation 0.5\n",
      "[0.10329022796802324, 0.02662299818877489, 0.007383799942783786] pressure error\n",
      "[0.10417700714847174, 0.049037576768854035, 0.016113648245047286] flux error\n",
      "Triangular errors - no perturbations\n",
      "[0.04865380618959045, 0.012528781793590854, 0.0031466771286353342] pressure error\n",
      "[0.06787911458237876, 0.020085847188375933, 0.005705542214407857] flux error\n",
      "Triangular errors - perturbation 0.5\n",
      "[0.05695409168995138, 0.015068959494986567, 0.003871569621335528] pressure error\n",
      "[0.053652948523950576, 0.031155048071189857, 0.011882770637109161] flux error\n"
     ]
    }
   ],
   "source": [
    "# Analytical solution\n",
    "x, y = sympy.symbols('x y')\n",
    "u = x * (1-x) * sympy.sin(x) * sympy.cos(y)\n",
    "u_f = sympy.lambdify((x, y), u, 'numpy')\n",
    "dux = sympy.diff(u, x)\n",
    "duy = sympy.diff(u, y)\n",
    "dux_f = sympy.lambdify((x, y), dux, 'numpy')\n",
    "duy_f = sympy.lambdify((x, y), duy, 'numpy')\n",
    "rhs = -sympy.diff(dux, x) - sympy.diff(duy, y)\n",
    "rhs_f = sympy.lambdify((x, y), rhs, 'numpy')\n",
    "\n",
    "###\n",
    "# This is where parameters can be modified to alter the convergence test.\n",
    "# The remaining lines should \n",
    "np.random.seed(42)\n",
    "base = 4\n",
    "domain = np.array([1, 1])\n",
    "basedim = np.array([base, base])\n",
    "num_refs = 3\n",
    "grid_type = 'cart'\n",
    "\n",
    "def run_convergence(grid_type, pert):\n",
    "    u_err = []\n",
    "    flux_err = []\n",
    "\n",
    "    for g in setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain):\n",
    "        # Reset the random seed for every grid realization.\n",
    "        # This should make no difference for the convergence test, \n",
    "        # but it makes sure that we can run unit tests based on the values obtained\n",
    "        # here.\n",
    "        np.random.seed(42)\n",
    "        \n",
    "        # Permeability tensor\n",
    "        k = perm(2, np.ones(g.num_cells))\n",
    "\n",
    "        # Set type of boundary conditions - Dirichlet\n",
    "        bound_faces = g.tags['domain_boundary_faces'].nonzero()[0]\n",
    "        bound_cond = bc.BoundaryCondition(g, bound_faces, ['dir'] * bound_faces.size)\n",
    "        \n",
    "        # MPFA discretization, and system matrix\n",
    "        flux, bound_flux, _, _ = mpfa.mpfa(g, k, bound_cond)\n",
    "        div = fvutils.scalar_divergence(g)\n",
    "        a = div * flux\n",
    "        \n",
    "        # Boundary conditions\n",
    "        xf = g.face_centers\n",
    "        u_bound = np.zeros(g.num_faces)\n",
    "        u_bound[bound_faces] = u_f(xf[0, bound_faces], xf[1, bound_faces])\n",
    "        \n",
    "        # Right hand side - contribution from the solution and the boundary conditions\n",
    "        xc = g.cell_centers\n",
    "        b = rhs_f(xc[0], xc[1]) * g.cell_volumes - div * bound_flux * u_bound\n",
    "\n",
    "        # Solve system, derive fluxes\n",
    "        u_num = spsolve(a, b)\n",
    "        flux_num = flux * u_num + bound_flux * u_bound\n",
    "\n",
    "        # Exact solution\n",
    "        u_ex = u_f(xc[0], xc[1])\n",
    "        du_ex_faces = np.vstack((dux_f(xf[0], xf[1]), duy_f(xf[0], xf[1])))\n",
    "        flux_ex = -np.sum(g.face_normals[:2] * du_ex_faces, axis=0)\n",
    "        flux_diff = flux_num - flux_ex\n",
    "  \n",
    "        u_err.append(np.sqrt(np.sum(g.cell_volumes * (u_num - u_ex)**2)) /\n",
    "                     np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n",
    "        flux_err.append(np.sqrt(np.sum((g.face_areas ** g.dim) * flux_diff**2))/\n",
    "                        np.sqrt(np.sum((g.face_areas ** g.dim) * flux_ex**2)))\n",
    "    return u_err, flux_err\n",
    "\n",
    "grids = ['cart', 'triangular']\n",
    "\n",
    "\n",
    "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
    "print('Cartesian errors - no perturbations')\n",
    "print(u_cart_nopert,'pressure error')\n",
    "print(f_cart_nopert,'flux error')\n",
    "\n",
    "u_cart_pert, f_cart_pert = run_convergence('cart', 0.5)\n",
    "print('Cartesian errors - perturbation 0.5')\n",
    "print(u_cart_pert,'pressure error')\n",
    "print(f_cart_pert,'flux error')\n",
    "\n",
    "\n",
    "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n",
    "print('Triangular errors - no perturbations')\n",
    "print(u_triang_nopert,'pressure error')\n",
    "print(f_triang_nopert,'flux error')\n",
    "\n",
    "\n",
    "u_triang_pert, f_triang_pert = run_convergence('triangular', 0.5)\n",
    "print('Triangular errors - perturbation 0.5')\n",
    "print(u_triang_pert,'pressure error')\n",
    "print(f_triang_pert,'flux error')\n",
    "\n",
    "# EK: These values were hard-coded 11.06.2017.\n",
    "assert np.abs(u_cart_nopert[2] - 0.0065122960444455652) < 1e-10\n",
    "assert np.abs(f_cart_nopert[2] - 0.011371514390600204) < 1e-10\n",
    "assert np.abs(u_cart_pert[2] - 0.0073837999427839654) < 1e-10\n",
    "assert np.abs(f_cart_pert[2] - 0.016113648245047234) < 1e-10\n",
    "\n",
    "# RB: These values were hard-coded 06.09.2018.\n",
    "assert np.abs(u_triang_nopert[2] - 0.0031466771286353342) < 1e-10\n",
    "assert np.abs(f_triang_nopert[2] - 0.005705542214407857) < 1e-10\n",
    "assert np.abs(u_triang_pert[2] - 0.003871569621335528) < 1e-10\n",
    "assert np.abs(f_triang_pert[2] - 0.011882770637109161) < 1e-10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neumann boundary condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian errors - no perturbations\n",
      "[0.2177630338617197, 0.053854028644456195, 0.013418192291060563] pressure error\n",
      "[4.04357927 4.01350849] convergence rate u^n-1 / u^n\n",
      "\n",
      "[0.020449627547754923, 0.00558959841388443, 0.0014594922987497814] flux error\n",
      "[3.65851463 3.82982385] convergence rate f^n-1 / f^n\n",
      "\n",
      "Triangular errors - no perturbations\n",
      "[0.1154963072658415, 0.028265747189434955, 0.007020068836548777] pressure error\n",
      "[4.08608718 4.02642023] convergence rate u^n-1 / u^n\n",
      "\n",
      "[0.04223113257437278, 0.011882184577715638, 0.0033976662367782447] flux error\n",
      "[3.55415558 3.49716062] convergence rate f^n-1 / f^n\n"
     ]
    }
   ],
   "source": [
    "# Analytical solution\n",
    "x, y = sympy.symbols('x y')\n",
    "u = x * (1-x) * sympy.sin(x) * sympy.cos(y)\n",
    "u_f = sympy.lambdify((x, y), u, 'numpy')\n",
    "dux = sympy.diff(u, x)\n",
    "duy = sympy.diff(u, y)\n",
    "dux_f = sympy.lambdify((x, y), dux, 'numpy')\n",
    "duy_f = sympy.lambdify((x, y), duy, 'numpy')\n",
    "rhs = -sympy.diff(dux, x) - sympy.diff(duy, y)\n",
    "rhs_f = sympy.lambdify((x, y), rhs, 'numpy')\n",
    "\n",
    "###\n",
    "# This is where parameters can be modified to alter the convergence test.\n",
    "# The remaining lines should \n",
    "np.random.seed(42)\n",
    "base = 4\n",
    "domain = np.array([1, 1])\n",
    "basedim = np.array([base, base])\n",
    "num_refs = 3\n",
    "grid_type = 'cart'\n",
    "\n",
    "def run_convergence(grid_type, pert):\n",
    "    u_err = []\n",
    "    flux_err = []\n",
    "\n",
    "    for g in setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain):\n",
    "        # Reset the random seed for every grid realization.\n",
    "        # This should make no difference for the convergence test, \n",
    "        # but it makes sure that we can run unit tests based on the values obtained\n",
    "        # here.\n",
    "        np.random.seed(42)\n",
    "        \n",
    "        # Permeability tensor\n",
    "        k = perm(2, np.ones(g.num_cells))\n",
    "        \n",
    "        # Exact solution\n",
    "        xf = g.face_centers\n",
    "        xc = g.cell_centers\n",
    "        u_ex = u_f(xc[0], xc[1])\n",
    "        du_ex_faces = np.vstack((dux_f(xf[0], xf[1]), duy_f(xf[0], xf[1])))\n",
    "        flux_ex = -np.sum(g.face_normals[:2] * du_ex_faces, axis=0)\n",
    "\n",
    "        # Set type of boundary conditions - Dirichlet\n",
    "        bound_faces = g.tags['domain_boundary_faces'].nonzero()[0]\n",
    "        n = g.nodes\n",
    "        top = bc.face_on_side(g, 'ymax')\n",
    "        bot = bc.face_on_side(g, 'ymin')\n",
    "        left = bc.face_on_side(g, 'xmin')\n",
    "        right = bc.face_on_side(g, 'xmax')\n",
    "\n",
    "        dir_faces = np.asarray(right[0])\n",
    "        neu_faces = np.hstack((left[0], bot[0], top[0]))\n",
    "        bound_cond = bc.BoundaryCondition(g, dir_faces, ['dir'] * len(dir_faces))\n",
    "                \n",
    "        # MPFA discretization, and system matrix\n",
    "        flux, bound_flux, _, _ = mpfa.mpfa(g, k, bound_cond)\n",
    "        div = fvutils.scalar_divergence(g)\n",
    "        a = div * flux\n",
    "        \n",
    "        # Boundary conditions\n",
    "        nfi, _, sgn = sps.find(g.cell_faces[neu_faces,:])\n",
    "        u_bound = np.zeros(g.num_faces)\n",
    "        u_bound[dir_faces] = u_f(xf[0, dir_faces], xf[1, dir_faces])\n",
    "        # innflow is always positive, so need to flip flux according to normal direction.\n",
    "        u_bound[neu_faces[nfi]] = flux_ex[neu_faces[nfi]] * (sgn)\n",
    "\n",
    "        # Right hand side - contribution from the solution and the boundary conditions\n",
    "        xc = g.cell_centers\n",
    "        b = rhs_f(xc[0], xc[1]) * g.cell_volumes - div * bound_flux * u_bound\n",
    "\n",
    "        # Solve system, derive fluxes\n",
    "        u_num = spsolve(a, b)\n",
    "        flux_num = flux * u_num + bound_flux * u_bound\n",
    "\n",
    "        # Exact solution\n",
    "        flux_diff = flux_num - flux_ex\n",
    "  \n",
    "        u_err.append(np.sqrt(np.sum(g.cell_volumes * (u_num - u_ex)**2)) /\n",
    "                     np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n",
    "        flux_err.append(np.sqrt(np.sum((g.face_areas ** g.dim) * flux_diff**2))/\n",
    "                        np.sqrt(np.sum((g.face_areas ** g.dim) * flux_ex**2)))\n",
    "    return u_err, flux_err\n",
    "\n",
    "grids = ['cart', 'triangular']\n",
    "\n",
    "## No tests on perturbed grids - I'm too lazy to find faces of the respective boundaries\n",
    "\n",
    "\n",
    "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
    "u_cart_rate = np.array(u_cart_nopert)\n",
    "u_cart_rate = u_cart_rate[:-1] / u_cart_rate[1:]\n",
    "f_cart_rate = np.array(f_cart_nopert)\n",
    "f_cart_rate = f_cart_rate[:-1] / f_cart_rate[1:]\n",
    "print('Cartesian errors - no perturbations')\n",
    "print(u_cart_nopert,'pressure error')\n",
    "print(u_cart_rate, 'convergence rate u^n-1 / u^n\\n')\n",
    "print(f_cart_nopert,'flux error')\n",
    "print(f_cart_rate, 'convergence rate f^n-1 / f^n\\n')\n",
    "\n",
    "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n",
    "u_triang_rate = np.array(u_triang_nopert)\n",
    "u_triang_rate = u_triang_rate[:-1] / u_triang_rate[1:]\n",
    "f_triang_rate = np.array(f_triang_nopert)\n",
    "f_triang_rate = f_triang_rate[:-1] / f_triang_rate[1:]\n",
    "\n",
    "print('Triangular errors - no perturbations')\n",
    "print(u_triang_nopert,'pressure error')\n",
    "print(u_triang_rate, 'convergence rate u^n-1 / u^n\\n')\n",
    "print(f_triang_nopert,'flux error')\n",
    "print(f_triang_rate, 'convergence rate f^n-1 / f^n')\n",
    "# EK: These values were hard-coded 11.06.2017.\n",
    "assert np.abs(u_cart_nopert[2] - 0.01341819229106214) < 1e-10\n",
    "assert np.abs(f_cart_nopert[2] - 0.0014594922987503009) < 1e-10\n",
    "\n",
    "# RB: These values were hard-coded 06.09.2018.\n",
    "assert np.abs(u_triang_nopert[2] - 0.007020068836548777) < 1e-10\n",
    "assert np.abs(f_triang_nopert[2] - 0.0033976662367782447) < 1e-10\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Robin boundary condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian errors - no perturbations\n",
      "[0.562515011909662, 0.14015471392164033, 0.034984342689359696] pressure error\n",
      "[4.01352902 4.00621258] convergence rate u^n-1 / u^n\n",
      "\n",
      "[0.22376060444435023, 0.061834410597352465, 0.01641158802403413] flux error\n",
      "[3.61870684 3.76772866] convergence rate f^n-1 / f^n\n",
      "\n",
      "Triangular errors - no perturbations\n",
      "[0.11192090940414492, 0.02516849817177693, 0.005995729167675974] pressure error\n",
      "[4.44686483 4.19773767] convergence rate u^n-1 / u^n\n",
      "\n",
      "[0.03311849158106908, 0.00866805360137708, 0.002734187490448117] flux error\n",
      "[3.82075297 3.17024843] convergence rate f^n-1 / f^n\n"
     ]
    }
   ],
   "source": [
    "# Analytical solution\n",
    "x, y = sympy.symbols('x y')\n",
    "u = x * (1-x) * sympy.sin(x) * sympy.cos(y)\n",
    "u_f = sympy.lambdify((x, y), u, 'numpy')\n",
    "dux = sympy.diff(u, x)\n",
    "duy = sympy.diff(u, y)\n",
    "dux_f = sympy.lambdify((x, y), dux, 'numpy')\n",
    "duy_f = sympy.lambdify((x, y), duy, 'numpy')\n",
    "rhs = -sympy.diff(dux, x) - sympy.diff(duy, y)\n",
    "rhs_f = sympy.lambdify((x, y), rhs, 'numpy')\n",
    "\n",
    "###\n",
    "# This is where parameters can be modified to alter the convergence test.\n",
    "# The remaining lines should \n",
    "np.random.seed(42)\n",
    "base = 4\n",
    "domain = np.array([1, 1])\n",
    "basedim = np.array([base, base])\n",
    "num_refs = 3\n",
    "grid_type = 'cart'\n",
    "\n",
    "def run_convergence(grid_type, pert):\n",
    "    u_err = []\n",
    "    flux_err = []\n",
    "\n",
    "    for g in setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain):\n",
    "        # Reset the random seed for every grid realization.\n",
    "        # This should make no difference for the convergence test, \n",
    "        # but it makes sure that we can run unit tests based on the values obtained\n",
    "        # here.\n",
    "        np.random.seed(42)\n",
    "        \n",
    "        # Permeability tensor\n",
    "        k = perm(2, np.ones(g.num_cells))\n",
    "        robin_weight = 1.2\n",
    "        # Exact solution\n",
    "        xf = g.face_centers\n",
    "        xc = g.cell_centers\n",
    "        u_ex = u_f(xc[0], xc[1])\n",
    "        du_ex_faces = np.vstack((dux_f(xf[0], xf[1]), duy_f(xf[0], xf[1])))\n",
    "        flux_ex = -np.sum(g.face_normals[:2] * du_ex_faces, axis=0)\n",
    "\n",
    "        # Set type of boundary conditions - Dirichlet\n",
    "        bound_faces = g.tags['domain_boundary_faces'].nonzero()[0]\n",
    "        n = g.nodes\n",
    "        top = bc.face_on_side(g, 'ymax')\n",
    "        bot = bc.face_on_side(g, 'ymin')\n",
    "        left = bc.face_on_side(g, 'xmin')\n",
    "        right = bc.face_on_side(g, 'xmax')\n",
    "\n",
    "        dir_faces = np.asarray(left[0])\n",
    "        rob_faces = np.hstack((right[0]))\n",
    "        neu_faces = np.hstack((bot[0], top[0]))\n",
    "        bc_faces = np.hstack([dir_faces, rob_faces])\n",
    "        bc_names = ['dir']*len(dir_faces) + ['rob'] * len(rob_faces)\n",
    "        bound_cond = bc.BoundaryCondition(g, bc_faces,  bc_names)\n",
    "                \n",
    "        # MPFA discretization, and system matrix\n",
    "        flux, bound_flux, _, _ = mpfa.mpfa(g, k, bound_cond, robin_weight=robin_weight)\n",
    "        div = fvutils.scalar_divergence(g)\n",
    "        a = div * flux\n",
    "        \n",
    "        # Boundary conditions\n",
    "        nfi, _, sgn_n = sps.find(g.cell_faces[neu_faces,:])\n",
    "        rfi, _, sgn = sps.find(g.cell_faces[rob_faces,:])\n",
    "        u_bound = np.zeros(g.num_faces)\n",
    "        u_bound[dir_faces] = u_f(xf[0, dir_faces], xf[1, dir_faces])\n",
    "        # innflow is always negative, so need to flip flux according to normal direction.\n",
    "        u_bound[neu_faces[nfi]] = flux_ex[neu_faces[nfi]] * (sgn_n)\n",
    "        u_bound[rob_faces[rfi]] = flux_ex[rob_faces[rfi]] * (sgn) +\\\n",
    "                                  robin_weight * u_f(xf[0, rob_faces], xf[1, rob_faces]) * g.face_areas[rob_faces]\n",
    "\n",
    "        # Right hand side - contribution from the solution and the boundary conditions\n",
    "        xc = g.cell_centers\n",
    "        b = rhs_f(xc[0], xc[1]) * g.cell_volumes - div * bound_flux * u_bound\n",
    "\n",
    "        # Solve system, derive fluxes\n",
    "        u_num = spsolve(a, b)\n",
    "        flux_num = flux * u_num + bound_flux * u_bound\n",
    "\n",
    "        # Exact solution\n",
    "        flux_diff = flux_num - flux_ex\n",
    "  \n",
    "        u_err.append(np.sqrt(np.sum(g.cell_volumes * (u_num - u_ex)**2)) /\n",
    "                     np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n",
    "        flux_err.append(np.sqrt(np.sum((g.face_areas ** g.dim) * flux_diff**2))/\n",
    "                        np.sqrt(np.sum((g.face_areas ** g.dim) * flux_ex**2)))\n",
    "    return u_err, flux_err\n",
    "\n",
    "grids = ['cart', 'triangular']\n",
    "\n",
    "## No tests on perturbed grids - I'm too lazy to find faces of the respective boundaries\n",
    "\n",
    "\n",
    "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
    "u_cart_rate = np.array(u_cart_nopert)\n",
    "u_cart_rate = u_cart_rate[:-1] / u_cart_rate[1:]\n",
    "f_cart_rate = np.array(f_cart_nopert)\n",
    "f_cart_rate = f_cart_rate[:-1] / f_cart_rate[1:]\n",
    "print('Cartesian errors - no perturbations')\n",
    "print(u_cart_nopert,'pressure error')\n",
    "print(u_cart_rate, 'convergence rate u^n-1 / u^n\\n')\n",
    "print(f_cart_nopert,'flux error')\n",
    "print(f_cart_rate, 'convergence rate f^n-1 / f^n\\n')\n",
    "\n",
    "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n",
    "u_triang_rate = np.array(u_triang_nopert)\n",
    "u_triang_rate = u_triang_rate[:-1] / u_triang_rate[1:]\n",
    "f_triang_rate = np.array(f_triang_nopert)\n",
    "f_triang_rate = f_triang_rate[:-1] / f_triang_rate[1:]\n",
    "\n",
    "print('Triangular errors - no perturbations')\n",
    "print(u_triang_nopert,'pressure error')\n",
    "print(u_triang_rate, 'convergence rate u^n-1 / u^n\\n')\n",
    "print(f_triang_nopert,'flux error')\n",
    "print(f_triang_rate, 'convergence rate f^n-1 / f^n')\n",
    "\n",
    "# RB: These values were hard-coded 11.09.2018.\n",
    "assert np.abs(u_cart_nopert[2] - 0.034984342689359696) < 1e-10\n",
    "assert np.abs(f_cart_nopert[2] - 0.01641158802403413) < 1e-10\n",
    "assert np.abs(u_triang_nopert[2] - 0.005995729167675974) < 1e-10\n",
    "assert np.abs(f_triang_nopert[2] - 0.002734187490448117) < 1e-10\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heterogeneous permeability field\n",
    "Jump in upper right corner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian errors - no perturbations\n",
      "[[0.23370055 0.23370055 0.23370055]\n",
      " [0.05302929 0.05302929 0.05302929]\n",
      " [0.01295075 0.01295075 0.01295075]] pressure error\n",
      "[[0.11072073 0.11072073 0.11072073]\n",
      " [0.02617215 0.02617215 0.02617215]\n",
      " [0.00645454 0.00645454 0.00645454]] flux error\n",
      "Cartesian errors - perturbation 0.5\n",
      "[[0.28144397 0.3116494  0.28115807]\n",
      " [0.08134885 0.0741203  0.06868011]\n",
      " [0.01484212 0.01673688 0.01728718]] pressure error\n",
      "[[0.18233422 0.18698603 0.20433553]\n",
      " [0.05353656 0.05186143 0.05246204]\n",
      " [0.01687373 0.01673338 0.01674571]] flux error\n",
      "Triangular errors - no perturbations\n",
      "[[0.18810597 0.19650285 0.19475267]\n",
      " [0.0490365  0.04950331 0.0494196 ]\n",
      " [0.01250799 0.01255227 0.01253359]] pressure error\n",
      "[[0.09610297 0.0942635  0.09702186]\n",
      " [0.02395548 0.02417054 0.02485514]\n",
      " [0.00653148 0.00662373 0.00680908]] flux error\n",
      "Triangular errors - perturbation 0.5\n",
      "[[0.25769202 0.30829373 0.30552852]\n",
      " [0.08177972 0.07242441 0.06817465]\n",
      " [0.01520087 0.01618674 0.01645769]] pressure error\n",
      "[[0.13567926 0.12937504 0.13327836]\n",
      " [0.04232624 0.04084399 0.04092574]\n",
      " [0.01302646 0.01309025 0.01333882]] flux error\n"
     ]
    }
   ],
   "source": [
    "# Analytical solution\n",
    "\n",
    "kappa_list = [1e-6, 1, 1e6]\n",
    "\n",
    "### End of parameter definitions\n",
    "\n",
    "# Permeability tensor, scalar for simplicity\n",
    "\n",
    "def chi(x, y):\n",
    "    return np.logical_and(np.greater(x, 0.5), np.greater(y, 0.5))\n",
    "\n",
    "# Analytical solution\n",
    "x, y = sympy.symbols('x y')\n",
    "u = sympy.sin(2 * x * pi) * sympy.sin(2 * y * pi) \n",
    "u_f = sympy.lambdify((x, y), u, 'numpy')\n",
    "dux = sympy.diff(u, x)\n",
    "duy = sympy.diff(u, y)\n",
    "dux_f = sympy.lambdify((x, y), dux, 'numpy')\n",
    "duy_f = sympy.lambdify((x, y), duy, 'numpy')\n",
    "rhs = -sympy.diff(dux, x) - sympy.diff(duy, y)\n",
    "rhs_f = sympy.lambdify((x, y), rhs, 'numpy')\n",
    "\n",
    "\n",
    "###\n",
    "# This is where parameters can be modified to alter the convergence test.\n",
    "# The remaining lines should \n",
    "np.random.seed(42)\n",
    "base = 4\n",
    "domain = np.array([1, 1])\n",
    "basedim = np.array([base, base])\n",
    "num_refs = 3\n",
    "grid_type = 'cart'\n",
    "\n",
    "def run_convergence(grid_type, pert):\n",
    "    u_err = np.zeros((num_refs, len(kappa_list)))\n",
    "    flux_err = np.copy(u_err)\n",
    "    for iter1, g in enumerate(setup_grids.grid_sequence_fixed_lines(basedim, num_refs, grid_type, pert, subdom_func=chi)):\n",
    "        # Reset the random seed for every grid realization.\n",
    "        # This should make no difference for the convergence test, \n",
    "        # but it makes sure that we can run unit tests based on the values obtained\n",
    "        # here.\n",
    "        np.random.seed(42)\n",
    "        \n",
    "        for iter2, kappa in enumerate(kappa_list):\n",
    "            char_func_cells = chi(g.cell_centers[0], g.cell_centers[1]) * 1.\n",
    "            perm_vec = (1 - char_func_cells) + kappa * char_func_cells\n",
    "            # Permeability tensor\n",
    "            k = perm(2, perm_vec)\n",
    "            bound_faces = g.tags['domain_boundary_faces'].nonzero()[0]\n",
    "            bound_cond = bc.BoundaryCondition(g, bound_faces, ['dir'] * bound_faces.size)\n",
    "            flux, bound_flux, _, _ = mpfa.mpfa(g, k, bound_cond)\n",
    "\n",
    "            xc = g.cell_centers\n",
    "            xf = g.face_centers\n",
    "            char_func_bound = chi(xf[0, bound_faces], xf[1, bound_faces]) * 1\n",
    "            \n",
    "            u_bound = np.zeros(g.num_faces)\n",
    "            u_bound[bound_faces] = u_f(xf[0, bound_faces], xf[1, bound_faces]) \\\n",
    "                                / ((1 - char_func_bound) + kappa * char_func_bound)\n",
    "\n",
    "            div = fvutils.scalar_divergence(g)\n",
    "            a = div * flux\n",
    "\n",
    "            b = rhs_f(xc[0], xc[1]) * g.cell_volumes - div * bound_flux * u_bound\n",
    "\n",
    "            u_num = spsolve(a, b)\n",
    "            u_ex = u_f(xc[0], xc[1]) / ((1 - char_func_cells) + kappa * char_func_cells)\n",
    "            \n",
    "            flux_num = flux * u_num + bound_flux * u_bound\n",
    "            du_ex_faces = np.vstack((dux_f(xf[0], xf[1]), duy_f(xf[0], xf[1])))\n",
    "            flux_ex = -np.sum(g.face_normals[:2] * du_ex_faces, axis=0)\n",
    "            flux_diff = flux_num - flux_ex\n",
    "\n",
    "            u_err[iter1, iter2] = (np.sqrt(np.sum(g.cell_volumes * (u_num - u_ex)**2)) /\n",
    "                         np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n",
    "            flux_err[iter1, iter2] = (np.sqrt(np.sum((g.face_areas ** g.dim) * flux_diff**2))/\n",
    "                            np.sqrt(np.sum((g.face_areas ** g.dim) * flux_ex**2)))\n",
    "\n",
    "    return u_err, flux_err\n",
    "\n",
    "grids = ['cart', 'triangular']\n",
    "\n",
    "## No tests on perturbed grids - I'm too lazy to find faces of the respective boundaries\n",
    "\n",
    "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
    "print('Cartesian errors - no perturbations')\n",
    "print(u_cart_nopert,'pressure error')\n",
    "print(f_cart_nopert,'flux error')\n",
    "\n",
    "u_cart_pert, f_cart_pert = run_convergence('cart', 0.5)\n",
    "print('Cartesian errors - perturbation 0.5')\n",
    "print(u_cart_pert,'pressure error')\n",
    "print(f_cart_pert,'flux error')\n",
    "\n",
    "\n",
    "u_triang_nopert, f_triang_nopert = run_convergence('triangular', 0)\n",
    "print('Triangular errors - no perturbations')\n",
    "print(u_triang_nopert,'pressure error')\n",
    "print(f_triang_nopert,'flux error')\n",
    "\n",
    "\n",
    "u_triang_pert, f_triang_pert = run_convergence('triangular', 0.5)\n",
    "print('Triangular errors - perturbation 0.5')\n",
    "print(u_triang_pert,'pressure error')\n",
    "print(f_triang_pert,'flux error')\n",
    "\n",
    "# EK: These values were hard-coded 11.06.2017.\n",
    "assert np.abs(u_cart_nopert[2, 0] - 0.01295075) < 1e-7\n",
    "assert np.abs(f_cart_nopert[2, 0] - 0.00645454) < 1e-7\n",
    "assert np.abs(u_cart_pert[2, 0] - 0.01484212) < 1e-7\n",
    "assert np.abs(f_cart_pert[2, 0] - 0.01687373) < 1e-7\n",
    "\n",
    "# RB: These values were hard-coded 06.09.2018.\n",
    "assert np.abs(u_triang_nopert[2, 0] - 0.01250799) < 1e-7\n",
    "assert np.abs(f_triang_nopert[2, 0] - 0.00653148) < 1e-7\n",
    "assert np.abs(u_triang_pert[2, 0] - 0.01520087) < 1e-7\n",
    "assert np.abs(f_triang_pert[2, 0] - 0.01302646) < 1e-7\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3d test cases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian errors - no perturbations\n",
      "[0.0038211823411524516, 0.0011419336515988462, 0.00030564723918128137] pressure error\n",
      "[0.03261861955556977, 0.010014070049802667, 0.002872556762934246] flux error\n",
      "Cartesian errors - perturbation 0.5\n",
      "[0.004190264016265817, 0.0012081169906913614, 0.0003273282405744348] pressure error\n",
      "[0.028005576371853302, 0.01068845032035332, 0.004161378449699483] flux error\n",
      "Tetrehedral errors - no perturbations\n",
      "[0.0019052595279009587, 0.00048473242101445527, 0.00012116396539578085] pressure error\n",
      "[0.01802429803935019, 0.007938868913250801, 0.0036882937583153238] flux error\n",
      "Tetrahedral errors - perturbation 0.3\n",
      "[0.002253246419480656, 0.0005353463651070588, 0.00013484254795268383] pressure error\n",
      "[0.019027005210865585, 0.009029346065771729, 0.004055069442646703] flux error\n"
     ]
    }
   ],
   "source": [
    "# Analytical solution\n",
    "x, y, z = sympy.symbols('x y z')\n",
    "u = sympy.cos(x) * sympy.sin(z) * sympy.cosh(y)\n",
    "u_f = sympy.lambdify((x, y, z), u, 'numpy')\n",
    "dux = sympy.diff(u, x)\n",
    "duy = sympy.diff(u, y)\n",
    "duz = sympy.diff(u, z)\n",
    "dux_f = sympy.lambdify((x, y, z), dux, 'numpy')\n",
    "duy_f = sympy.lambdify((x, y, z), duy, 'numpy')\n",
    "duz_f = sympy.lambdify((x, y, z), duz, 'numpy')\n",
    "rhs = -sympy.diff(dux, x) - sympy.diff(duy, y) - sympy.diff(duz, z)\n",
    "rhs_f = sympy.lambdify((x, y, z), rhs, 'numpy')\n",
    "\n",
    "###\n",
    "# This is where parameters can be modified to alter the convergence test.\n",
    "# The remaining lines should \n",
    "np.random.seed(42)\n",
    "base = 4\n",
    "domain = np.array([1, 1, 1])\n",
    "basedim = np.array([base, base, base])\n",
    "num_refs = 3\n",
    "\n",
    "def run_convergence(grid_type, pert):\n",
    "    u_err = []\n",
    "    flux_err = []\n",
    "\n",
    "    for g in setup_grids.grid_sequence(basedim, num_refs, grid_type, pert, domain=domain):\n",
    "        # Reset the random seed for every grid realization.\n",
    "        # This should make no difference for the convergence test, \n",
    "        # but it makes sure that we can run unit tests based on the values obtained\n",
    "        # here.\n",
    "        np.random.seed(42)\n",
    "        \n",
    "        # Permeability tensor\n",
    "        k = perm(3, np.ones(g.num_cells))\n",
    "\n",
    "        # Set type of boundary conditions - Dirichlet\n",
    "        bound_faces = g.tags['domain_boundary_faces'].nonzero()[0]\n",
    "        bound_cond = bc.BoundaryCondition(g, bound_faces, ['dir'] * bound_faces.size)\n",
    "        \n",
    "        # MPFA discretization, and system matrix\n",
    "        flux, bound_flux ,_,_= mpfa.mpfa(g, k, bound_cond)\n",
    "        div = fvutils.scalar_divergence(g)\n",
    "        a = div * flux\n",
    "\n",
    "        # Boundary conditions\n",
    "        xf = g.face_centers\n",
    "        u_bound = np.zeros(g.num_faces)\n",
    "        u_bound[bound_faces] = u_f(xf[0, bound_faces], xf[1, bound_faces], xf[2, bound_faces])\n",
    "        \n",
    "        # Right hand side - contribution from the solution and the boundary conditions\n",
    "        xc = g.cell_centers\n",
    "        b = rhs_f(xc[0], xc[1], xc[2]) * g.cell_volumes - div * bound_flux * u_bound\n",
    "\n",
    "        # Solve system, derive fluxes\n",
    "        u_num = spsolve(a, b)\n",
    "        flux_num = flux * u_num + bound_flux * u_bound\n",
    "\n",
    "        # Exact solution\n",
    "        u_ex = u_f(xc[0], xc[1], xc[2])\n",
    "        du_ex_faces = np.vstack((dux_f(xf[0], xf[1], xf[2]),\n",
    "                                 duy_f(xf[0], xf[1], xf[2]),\n",
    "                                 duz_f(xf[0], xf[1], xf[2])))\n",
    "        flux_ex = -np.sum(g.face_normals * du_ex_faces, axis=0)\n",
    "        flux_diff = flux_num - flux_ex\n",
    "\n",
    "        u_err.append(np.sqrt(np.sum(g.cell_volumes * (u_num - u_ex)**2)) /\n",
    "                     np.sqrt(np.sum(g.cell_volumes * u_ex**2)))\n",
    "        flux_err.append(np.sqrt(np.sum((g.face_areas ** g.dim) * flux_diff**2))/\n",
    "                        np.sqrt(np.sum((g.face_areas ** g.dim) * flux_ex**2)))\n",
    "    return u_err, flux_err\n",
    "\n",
    "u_cart_nopert, f_cart_nopert = run_convergence('cart', 0)\n",
    "print('Cartesian errors - no perturbations')\n",
    "print(u_cart_nopert,'pressure error')\n",
    "print(f_cart_nopert,'flux error')\n",
    "\n",
    "u_cart_pert, f_cart_pert = run_convergence('cart', 0.5)\n",
    "print('Cartesian errors - perturbation 0.5')\n",
    "print(u_cart_pert,'pressure error')\n",
    "print(f_cart_pert,'flux error')\n",
    "\n",
    "\n",
    "u_triang_nopert, f_triang_nopert = run_convergence('tetrahedral', 0)\n",
    "print('Tetrehedral errors - no perturbations')\n",
    "print(u_triang_nopert,'pressure error')\n",
    "print(f_triang_nopert,'flux error')\n",
    "\n",
    "u_triang_pert, f_triang_pert = run_convergence('tetrahedral', 0.3)\n",
    "print('Tetrahedral errors - perturbation 0.3')\n",
    "print(u_triang_pert,'pressure error')\n",
    "print(f_triang_pert,'flux error')\n",
    "\n",
    "# EK: These values were hard-coded 11.06.2017.\n",
    "assert np.abs(u_cart_nopert[2] -0.00030564723918128137 ) < 1e-10\n",
    "assert np.abs(f_cart_nopert[2] - 0.0028725567629342458) < 1e-10\n",
    "assert np.abs(u_cart_pert[2] - 0.00032732824057441688) < 1e-10\n",
    "assert np.abs(f_cart_pert[2] - 0.0041613784496996173) < 1e-10\n",
    "\n",
    "# RB: These values were hard-coded 06.09.2018.\n",
    "assert np.abs(u_triang_nopert[2] - 0.00012116396539578085) < 1e-10\n",
    "assert np.abs(f_triang_nopert[2] - 0.0036882937583153238) < 1e-10\n",
    "assert np.abs(u_triang_pert[2] - 0.00013484254795268383) < 1e-10\n",
    "assert np.abs(f_triang_pert[2] - 0.004055069442646703) < 1e-10"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
